Journal of Mathematical Biology

We introduce a spectral existence criterion for the evolution of cooperation in the form of the inequality{lambda} maxb > c, where{lambda} max is the leading eigenvalue of an interaction operator encoding population structure, and b and c represent benefit and cost tradeoffs, respectively. Nowaks five rules for the evolution of cooperation correspond to cases in which the cooperation condition reduces to a scalar assortment coefficient. These results follow from the Price equation, which sheds light on a long-standing debate on the role of inclusive fitness and evolutionary dynamics in explaining the evolution of cooperation.

Clavicle fractures often exhibit markedly different clinical outcomes: some patients recover acceptable function despite shortening or displacement, whereas others with apparently similar deformity develop persistent pain, functional loss, or poor healing. To explain this distinction, we propose a minimal nonlinear mechanical model for prognostic analysis of clavicle fractures. The model describes the interaction between fracture-related shortening and compensatory shoulder-girdle posture through a reduced equilibrium equation incorporating stiffness, geometric nonlinearity, and shortening-posture coupling. Within this framework, we analyze equilibrium branches, local stability, and the emergence of critical thresholds. We show that post-fracture destabilization can be interpreted as a fold bifurcation, while more complex parameter dependence gives rise to cusp-type structures and multistability. These bifurcation mechanisms provide a mathematical explanation for sudden deterioration after injury or treatment, as well as for strong inter-individual variability. We further introduce an optimization principle based on a utility functional to guide treatment planning. The analysis predicts that the optimal safe correction should lie strictly below the bifurcation threshold, thereby generating a natural safety margin. Although the model is simplified and has not yet been calibrated against patient data, it nevertheless provides a theoretical framework for understanding why fracture prognosis may deteriorate abruptly near critical mechanical conditions and offers a dynamical-systems interpretation of empirical treatment thresholds used in clinical practice.

The basic reproduction number of an infectious disease is known to depend on the structure of contacts between individuals in a population. This relationship has been explored mathematically through two well-known models: one which depends on a matrix of contact rates between different demographic groups, and another which depends on the variability of contact rates over the population. Here we introduce a model that combines and generalises these two approaches. We derive a formula for the basic reproduction number and validate it through comparisons to simulated outbreaks. Applying this method to contact survey data collected in Belgium between 2020 and 2022, we find that our model produces higher estimates of the basic reproduction number and larger relative changes over periods when social contact behaviour was changing during the COVID-19 pandemic. Our analysis suggests some practical considerations when using contact data in models of infectious disease transmission.

We introduce a mechanistic, nonlocal tumour-growth model designed specifically to capture explosive dynamics that are not adequately explained by standard logistic reaction-diffusion descriptions. The motivation is empirical: the universal scaling law reported in [1] provides compelling cross-sectional evidence of superlinear tumour activity versus tumour burden, but as a phenomenological relationship it does not by itself supply a dynamical mechanism, nor does it rigorously describe how explosive growth emerges, how fast it develops, or how spatial interactions and tissue boundaries influence it. Our model addresses this gap by incorporating nonlocal proliferative feedback--cells respond to a spatially aggregated neighbourhood signal--and a singular, Kawarada-type acceleration that produces "quenching": tumour density stays bounded while the proliferative drive becomes unbounded as the aggregated signal approaches a critical threshold. This offers a concrete mechanistic route to explosive escalation consistent with physical boundedness. We analyse the model under no-flux (Neumann) boundary conditions, appropriate for reflecting tissue interfaces. In the spatially homogeneous setting we prove finite-time onset of the explosive regime and obtain explicit rates for how rapidly it is approached. For spatially heterogeneous perturbations we derive a transparent spectral stability theory showing how the interaction kernel selects spatial scales and how the singular acceleration tightens stability margins as the explosive threshold is approached. These results provide interpretable links between nonlocal interaction structure, boundary effects, and the emergence of rapid growth. Finally, to connect mechanism to data in the spirit of [1], we embed the model in a Bayesian inference framework that treats the interaction kernel and the acceleration strength as unknown and learned from tumour-growth observations. This enables uncertainty-aware estimation of explosive onset times, escalation rates, and stability margins, while positioning the scaling law of [1] as an observable signature that our mechanistic model can explain and quantify rather than merely fit.

Cancer-associated fibroblasts (CAFs) are a central component of the tumor microenvironment that facilitate a supportive niche for cancer progression and metastasis. Experimental evidence suggests that CAFs can facilitate estrogen-independent tumor growth, thereby reducing the efficacy of anti-hormonal therapies. Understanding and quantifying the complex interactions between tumor cells, hormonal signalling, and the microenvironment are crucial for designing more effective and individualized treatment strategies. We propose a mathematical framework to explore the influence of CAFs on ER+ breast cancer progression and to evaluate strategies to mitigate their impact. We develop a system of nonlinear ordinary differential equations that substantiates the experimental observations by providing a mechanistic basis for the role of CAFs in regulating estrogen-independent growth dynamics. We then employ optimal control theory to evaluate distinct therapeutic approaches involving monotherapy or combinations of: (i) inhibition of tumor-to-CAF signaling, (ii) inhibition of CAF-to-tumor proliferative signaling, and (iii) endocrine therapy. Taken together, our results demonstrate that CAF-targeted strategies can enhance treatment efficacy across various estrogen dosing regimens. Our study provides new insights into the potential of CAF as a therapeutic target that could help to improve existing approaches for endocrine therapies.

Cytokine-mediated communication is a central mechanism by which immune cells coordinate activation, differentiation and proliferation. While mechanistic reaction-diffusion models provide detailed descriptions of cytokine secretion and uptake at the cellular scale, their computational cost limits their applicability to large and densely packed cell populations. Previously employed approximations of cytokine diffusion fields rely on assumptions that neglect the influence of cellular geometry and volume exclusion. In this work, we study a macroscopic description of cytokine diffusion and reaction dynamics based on homogenization techniques, rigorously linking microscopic reaction-diffusion formulations to effective continuum models. The resulting homogenized equations replace discrete responder cells with a continuous density, while retaining essential features of cellular uptake and excluded-volume effects. Further, we show that in regimes with approximate radial symmetry, classical Yukawa-type solutions emerge as limiting cases of the homogenized model, provided appropriate correction factors are included. Overall, our approach allows efficient multiscale modeling of cytokine signaling in complex immune-cell environments.

BackgroundTo control leishmaniasis, chemotherapy drugs are currently under development. However, these drugs often exhibit poor efficacy and are associated with toxicity, adverse effects, and drug resistance. At present, no specific drug is available for the treatment of leishmaniasis. Meanwhile, vaccine research is ongoing. Recent studies have analysed some experimental vaccines using mathematical models. AimIn previous work, drug targeting was focused on the entire human body rather than specifically addressing infected macrophages and parasites. In our current approach, we aim to eliminate infected macrophages and parasites through nano-drug design. Specifically, we utilise two types of nanoparticles: iron oxide and citric acid-coated iron oxide. Moving forward, we plan to advance this strategy using mathematical modelling of macrophage-parasite interactions. MethodsWe design PDE-based models of macrophages and parasites, incorporating cytokine dynamics, to support nano-drug development. Drug efficacy is estimated using posterior distributions to analyse phenotypic fluctuations of macrophages and parasites during the design phase. We investigate implicit and semi-implicit treatment schemes, focusing on energy decay properties. To model drug flow during treatment, we introduce a three-phase moving boundary problem. Comparative analyses are conducted to evaluate macrophage and parasite behaviour with and without treatment. Finally, the entire framework is implemented within a virtual lab environment. ResultsThe results show that the nano-drug exhibits better efficacy compared to combined drug doses. We analysed and compared two types of nano-drug particles: iron oxide and citric acid-coated iron oxide. We discuss how the drug effectively targets and eliminates infected macrophages and parasites. ConclusionOur models results and simulations will support researchers conducting further studies in nano-drug design for leishmaniasis. These simulations are performed within a virtual lab environment.

The environments inhabited by flying insects demand a balance between flight efficiency and flight manoeuvrability. In structural oscillators such as the insect indirect flight motor, efficiency (arising from resonance) and manoeuvrability (arising from kinematic modulation) are typically quid pro quo, with modulation incurring penalties to efficiency. Band-type resonance is a phenomenon that offers, in theory, a strategy to lessen these penalties via careful navigation through a band of efficient kinematic states. However, identifying this band is challenging: no methods exist to identify the complete band in realistic motor models, involving elasticity distributed across thorax and wing. Nor are the effects of elasticity distribution on the band known. In this work, we address both open topics. We present a suite of numerical methods for identifying the complete resonance band in general systems. Applying them to models of the insect flight motor with distributed elasticity--thoracic and wing flexion--reveals that distributed elasticity is moderate-risk but high-reward morphological feature. Well-tuned distributions expand the resonance band over fourfold whereas poorly-tuned distributions completely extinguish the resonance band. These results indicate that distributing elasticity across the insect flight motor can have adaptive value, and motivate broader work identifying distributions across species.

Mitochondria often form branching membrane networks distributed throughout the cell interior. In many, though not all, cell types, these networks are observed to consist of one large connected component together with many smaller fragments. Why does this pattern arise? Does it reflect a specific biological function, an external biophysical constraint, or something simpler? Using results from extremal graph theory, we prove a new theorem which suggests that, under a sufficiently broad sampling of the space of mitochondria-like graphs, the predominance of three-way junctions makes the appearance of a large component likely. This suggests that, in some settings, a large component may serve as a useful null model for mitochondrial network structure rather than requiring a dedicated explanation. More broadly, our result points towards testable predictions, since systematic deviations from this baseline may help reveal additional constraints or mechanisms shaping mitochondrial morphology.

We describe an approach for analyzing biological networks using rows of the Krylov subspace of the adjacency matrix. Specifically, we explore the scenario where the Krylov subspace matrix is computed via power iteration using a non-random and potentially non-uniform initial vector that captures a specific biological state or perturbation. In this case, the rows the Krylov subspace matrix (i.e., Krylov trajectories) carry important functional information about the network nodes in the biological context represented by the initial vector. We demonstrate the utility of this approach for community detection and perturbation analysis using the C. Elegans neural network.

Recent progress in mathematical kinship modelling has allowed one to predict the probable numbers of kin for a typical population member. In the models, kin may be structured by age and sex, both in static or time-variant demographies. Knowing the probable numbers of kin in different stages - such as parity, health status, or geographic location - however, remains an open challenge in Kinship Demography. Knowing how population structure delimits kin to distinct stages is an advance - for instance, the probability of having one sister at home and one sister away has different social implications from the probability of having two sisters. We present a novel analytical framework, grounded in branching process theory, that provides kin-number distributions jointly structured by age and stage. Using recursive compositions of probability generating functions (PGFs), we derive the joint age, stage, and age x stage kin-number distributions. All marginal distributions over either dimension naturally emerge. Simple extensions of the PGF approach additionally yield: the joint distribution of an individuals own stage and their kins stage; the probable numbers of kin deaths, both in total and by generation number; and the probabilities of being kinless and/or orphaned. We demonstrate the framework through novel results in an application using UK parity-specific fertility and mortality data. HighlightsO_LIA new method calculates probability generating functions for the number of kin structured by age and stage C_LIO_LIThe model allows predicting the probable numbers of kin organised by age and stage C_LIO_LIRecursive nesting of probability generating functions in branching processes is used C_LIO_LIAn application is presented highlighting the novel results C_LI

How genetically identical cells spontaneously break symmetry to assume divergent fates is a fundamental problem in developmental biology. While modern genomics has mapped the vast molecular repertoire involved in gene regulation, understanding the mechanism of cell state transitions that drive differentiation remains a formidable challenge. To address this, we use a reaction-kinetic framework to analyze recurring motifs of two and three competing master regulators. While typically such circuits are studied numerically, we show that assuming symmetry in nodes and interactions provides exact analytical description of the bifurcations governing cell fate transitions. We find that the possible cell fates across all considered topologies are dictated by a single dimensionless quantity, {beta}--the ratio of protein degradation to production rates. In the binary Toggle Switch (TS), decreasing {beta} destabilizes the symmetric (stem cell) state, giving rise to two asymmetric (differentiated) fates via a supercritical pitchfork bifurcation. In the three-component Toggle Triad (TT), low values of {beta} yield three asymmetric fates through subcritical pitchfork bifurcation, creating an intermediate range of {beta} where both symmetric and asymmetric fates are simultaneously stable. For the Self-Activating Toggle Switch (SATS), we identify a new parameter for the self-activation threshold ({theta}) and show that decreasing{theta} progressively stabilizes the uncommitted state, leading to a regime of tristability. Building on these temporal bifurcations, we next address the feasibility of spatial structure formation: can these multistable fates stably coexist within a spatial domain? Through a minimal model of cell-cell communication via free diffusion, we extend these motifs into reaction-diffusion systems, which reveals a direct role of network topology on spatial organization. We prove that any heterogeneous pattern in two-node circuits is inherently transient and unstable. In contrast, the three-node repressive network supports the stable spatial coexistence of differentiated phenotypes through pure diffusion, a phenomenon we analyze by studying heteroclinic interface solutions as building blocks. By reducing complex regulatory dynamics to tractable models with physically meaningful parameters, we establish a minimal framework which relates topology to cell fate. Finally, the effects of temporal multistability on pattern formation provide an excellent studying ground for morphogenesis, synthetic biology, and the overarching problem of spatiotemporal self-organization.

Human walking is intrinsically variable. For example, there is considerable stride to stride variability even when walking speed is constant. This variability is due to uncertainty in the sensorimotor system and the environment, and is shaped by both musculoskeletal dynamics (e.g. joint stiffness and damping originating from muscles) and the control strategy used to mitigate the effects of uncertainty. Yet, insight into how sensorimotor noise shapes walking variability is limited due to a lack of experimental methods to assess sensorimotor noise and control strategies during walking. Simulations that account for uncertainty can elucidate how sensorimotor noise affects movement variability but due to numerical challenges, accounting for sensorimotor noise is not common in simulations of walking. Existing simulations have hugely simplified musculoskeletal dynamics (e.g. no muscles), the control policy (e.g. pre-defined feedback loops), or sensorimotor noise sources (e.g. only motor noise). Here, we performed stochastic optimal control simulations of walking based on a model with 9 degrees of freedom and 18 muscles to study how the level of sensory and motor noise influences walking. We solved for feedforward muscle excitations and full-state time-varying feedback gains that minimised expected effort while generating periodic, and hence stable, gait patterns. To enable these simulations, we approximated the state distribution with a Gaussian and used an unscented transform to propagate the state covariance. Resulting optimisation problems were solved with direct collocation. Sensorimotor noise level had a small effect on the mean kinematics but shaped kinematic and muscle activity variability as well as expected effort. Although simulations underestimated the magnitude of experimental positional variability, they captured its structure. In agreement with experimental results, the control policy prioritised limiting variability of centre of mass kinematics and minimal swing foot clearance over limiting joint angle variability. Hence, our simulations suggest that effort minimisation underlies these observations. Author summaryWhen performing a movement multiple times, each repetition will be slightly different due to random disturbances in the neural signals used to control movement, i.e. sensorimotor noise. Because it is difficult to measure inside the nervous system of a moving person, computer simulations are used to study movement control. They found that both sensorimotor noise and musculoskeletal mechanics determine how people control arm movements and standing. However, there are no simulations of walking that systematically evaluated how sensorimotor noise level influences walking kinematics because they pose computational challenges. Here, we proposed and used an approach for minimal effort simulations of walking in the presence of uncertainty. We imposed forward speed and stability but not kinematics. We found that the level of sensorimotor noise had little effect on the mean movement but a strong effect on the variability and the expected effort. The control strategy prioritised reducing the variability of the centre of mass position and swing foot clearance over reducing the variability of individual joint angles, which is also observed in experiments. Interestingly, strict control of centre of mass position and foot clearance in our simulations emerged from minimising effort.

A key component of intraneuronal communication is the modulation of postsynaptic firing frequencies by stochastic transmitter release from presynaptic neurons. The time interval between successive postsynaptic firings is called the inter-spike interval (ISI), and understanding its statistics is integral to neural information processing. We start with a model of an excitatory chemical synapse with postsynaptic neuron firing governed as per a classical integrate-and-fire model. Using a first-passage time framework, we derive exact analytical results for the ISI statistical moments, revealing parameter regimes driving precision in postsynaptic action potential timing. Next, we extended this analysis to include both an excitatory and an inhibitory presynaptic connection onto the same postsynaptic neuron. We consider both a fixed postsynaptic-firing threshold and a threshold that adapts based on the postsynaptic membrane potential history. Our analysis shows that the latter adaptive threshold can result in scenarios where increasing the inhibitory input frequency increases the postsynaptic firing frequency. Moreover, we characterize parameter regimes where ISI noise is hypo-exponential or hyperexponential based on its coefficient of variation being less than or higher than one, respectively.

Mass-conserved reaction-diffusion systems are used as mathematical models for various phenomena such as cell polarity. Numerical simulations of this system present transient dynamics in which multiple stripe patterns converge to spatially monotonic patterns. Previous studies indicated that the transient dynamics are driven by a mass conservation law and by variations in the amount of substance contained in each pattern, which we refer to as "pattern flux". However, it is challenging to mathematically investigate these pattern dynamics. In this study, we introduce a reaction-diffusion compartment model to investigate the pattern dynamics in view of the conservation law and the pattern flux. This model is defined on multiple intervals (compartments), and diffusive couplings are imposed on each boundary of the compartments. Corresponding to the transient dynamics in the original system, we consider the dynamics around stripe patterns in the compartment model. We derive ordinary differential equations describing the pattern dynamics of the compartment model and analyze the existence and stability of equilibria for the reduced ODE with respect to the boundary parameters. For a specific parameter setting, we obtained results consistent with previous studies. Moreover, we present that the stripe patterns in the compartment model are potentially stabilized by changing the parameter, which is not observed in the original system. We expect that the methodology developed in this paper is extendable to various directions, such as membrane-induced pattern control.

Lightweight epidemic calculators are widely used for teaching and rapid scenario exploration, yet many omit the methodological detail needed for scientific reuse. We present a browser-native SIR calculator that exposes forward Euler and classical fourth-order Runge--Kutta (RK4) integration alongside epidemiologically interpretable outputs and a population-conservation diagnostic. The implementation is anchored to analytical properties of the deterministic SIR system, including the epidemic threshold, the peak condition, and the final-size relation. Benchmark experiments show that RK4 is essentially step-size invariant over practical discretizations, whereas Euler at a coarse one-day step overestimates peak prevalence by 3.97% and final size by 0.66% relative to a fine-step RK4 reference. These results demonstrate that browser-based tools can support publication-quality computational narratives when solver choice, diagnostics, and assumptions are treated as first-class outputs.

Neuronal activity alters extracellular ion concentrations and electric potentials. Ephaptic effects refer to the feedback influence that these extracellular changes can have on neuronal activity. While electric ephaptic effects occur on a fast timescale due to extracellular potential perturbations, ionic ephaptic effects are driven by slower, accumulative changes in ion concentrations. Among the previous computational studies of ephaptic effects, the vast majority have focused exclusively on electric effects, while ionic ephaptic effects have largely been neglected. In this work, we present an electrodiffusive computational framework consisting of two-compartment neurons that interact via a shared extracellular space. By accounting for both electric potentials and ion-concentration dynamics in a self-consistent manner, our framework enables us to explore the relative roles of electric and ionic ephaptic effects. Through numerical experiments, we demonstrate that ionic and electric ephaptic interactions play very different roles. While ionic ephaptic interactions increase population firing rates, electric ephaptic interactions primarily drive subtle shifts in spike timing. Furthermore, we show that these spike shifts cause the phase difference (the distance in spike times between a small collection of neurons) to converge to a stable, unique phase difference, which we coin the ephaptic intrinsic phase preference. Author summaryNeurons predominantly communicate through synapses: specialized contact points where a brief electrical signal, known as a spike or action potential, in one neuron influences another. Neurons generate these spikes by exchanging ions with the surrounding extracellular space. This way, spiking neurons alter extracellular ion concentrations and electric potentials. Since neurons are sensitive to such changes in their environment, they can also influence one another indirectly through the shared extracellular medium. This form of non-synaptic interaction is known as ephaptic coupling. Most computational models of neuronal activity neglect ephaptic interactions, and those that include them typically consider only electric effects while ignoring ionic contributions. As a result, the relative roles of electric and ionic ephaptic effects remain poorly understood. Here, we introduce a computational framework that accounts for both mechanisms in a self-consistent way. Our results show a functional distinction: ionic ephaptic effects act slowly, regulating population firing rates, whereas electric ephaptic effects act on millisecond timescales and subtly shift spike timing. These shifts cause spike-time differences between neurons to converge to a stable value, a phenomenon we call ephaptic intrinsic phase preference.

Evolving populations, especially in the strong-selection-weak-mutation limit, can be modelled as adaptive walks on fitness landscapes, moving in fitness-increasing mutational steps until reaching a fitness peak--a local optimum. Simulations of such adaptive walks--on a multi-peaked empirical landscape of the folA gene and on landscapes generated by the Rough Mount Fuji (RMF) model-- have shown that some landscapes are highly navigable, meaning that the highest x% of peaks are reached by >> x% of adaptive walks. This prompts the question of how adaptive walks can be so successful despite the local, myopic rules behind each adaptive step. Here, we investigate this question using simulations and mathematical approximations of random adaptive walks on a simplified RMF landscape. The landscape has a low-to-intermediate fitness region, whose size reconciles a low peak density with a high peak number. Despite the high number of peaks, walkers are likely to exit this region without terminating at a peak because the probability of a peak transition at each step is low and a fitness gradient guides walkers to the high-fitness region in few steps. Thus, three features are sufficient to explain why adaptive walks in the simplified RMF landscape are likely to reach a small fraction of top-ranking peaks: a low-to-intermediate fitness region with a high number of peaks, a low peak-transition probability, and which is crossed in few steps. We find that these three features are also present in the empirical folA landscape, suggesting that similar principles may apply.

At the fundamental conceptual level, two alternatives have traditionally been considered for how mutations arise and how evolution happens: 1) random mutation and natural selection, and 2) Lamarckism. Recently, the theory of Interaction-based Evolution (IBE) has been proposed, according to which mutations are neither random nor Lamarckian, but are influenced by information accumulating internally in the genome over generations. Based on the estimation-of-distribution algorithms framework, we present a simulation model that demonstrates nonrandom, non-Lamarckian mutation concretely while capturing indirectly several aspects of IBE: selection, recombination, and nonrandom, non-Lamarckian mutation interact in a complementary fashion; evolution is driven by the interaction of parsimony and fit; and random bits do not directly encode improvement but enable generalization by the manner in which they connect with the rest of the evolutionary process. Connections are drawn to Darwins observations that changed conditions increase the rate of production of heritable variation; to the causes of bell-shaped distributions of traits and how these distributions respond to selection; and to computational learning theory, where analogizing evolution to learning in accord with IBE casts individuals as examples and places the learned hypothesis at the population level. The model highlights the importance of incorporating internal integration of information through heritable change in both evolutionary theory and evolutionary computation.

Agent-based modeling (ABM) is a computational method for predicting the emergent outcomes of interacting, autonomous individuals in a complex system. Here, ABM is used to simulate interactions between fibroblast and myofibroblast cells during idiopathic pulmonary fibrosis (IPF) in alveolar tissue microenvironments. These microenvironments are derived from histology of a healthy human lung sample and moderate- and severe-IPF lung samples. Fibroblast differentiation, cell migration, and collagen secretion in response to the spatial distribution of the cytokine transforming growth factor-beta are captured in the ABM using NetLogo software. Results are presented from one simulated year without treatment and with mechanisms representing treatment by pirfenidone and pentoxifylline, alone and in combination. A total of 180 in silico experiments are run, analyzed, and compared in a high-throughput workflow. The effects of the initial number of fibroblasts and treatment scenarios on various metrics related to collagen accumulation and collagen invasion into alveolar regions are determined. The ABM and the analysis files are shared to facilitate model reuse. By integrating computational modeling of IPF and therapeutics, this research aims to improve understanding of fibrosis progression and assess the efficacy of novel and existing treatments targeting different mechanisms to inform decision-making for IPF treatment.